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We develop a new formalism to study nonlinear evolution in the growth of large-scale structure, 
by following the dynamics of gravitational clustering as it builds up in time. This approach is 
conveniently represented by Feynman diagrams constructed in terms of three objects: the initial 
5^ ■ conditions (e.g. perturbation spectrum), the vertex (describing non-linearities) and the propagator 

(describing linear evolution). We show that loop corrections to the linear power spectrum organize 
themselves into two classes of diagrams: one corresponding to mode-coupling effects, the other to a 
renormalization of the propagator. Resummation of the latter gives rise to a quantity that measures 
the memory of perturbations to initial conditions as a function of scale. As a result of this, we 
show that a well-defined (renormalized) perturbation theory follows, in the sense that each term 
in the remaining mode-coupling series dominates at some characteristic scale and is subdominant 
otherwise. This is unlike standard perturbation theory, where different loop corrections can become 
of the same magnitude in the nonlinear regime. In companion papers we compare the resummation 
of the propagator with numerical simulations, and apply these results to the calculation of the 
nonlinear power spectrum. Remarkably, the expressions in renormalized perturbation theory can 
be written in a way that closely resembles the halo model. 



I. INTRODUCTION: IMPROVING PERTURBATION THEORY 
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The growth of density and velocity perturbations is a fundamental result of cosmology that can be used 
to test basic properties of the universe, such as the amount of dark matter and dark energy. Although the 
linear and weakly non-linear growth of perturbations is well understood, giving the large-scale asymptotics 
(tree-level amplitudes in diagrammatic perturbation theory language [l|) of the power spectrum and higher- 
order statistics as a function of time, the nonlinear corrections to these results ("loop corrections" 0) are 
less well understood (see 3] for a review). 

The standard approach of perturbation theory (hereafter PT) is to solve the equations of motion by 
expanding in powers of the perturbation amplitudes. This is well justified when doing tree-level PT, as 
one is essentially looking for asymptotic behavior (the large-scale limit, where perturbation amplitudes 
become very small compared to unity), however as one approaches the nonlinear scale and loop corrections 
become increasingly important, the validity of such a method becomes less clear. Part of the issue is that 
cosmological PT is not a standard PT in the sense of having a small coupling constant; its validity depends 
on the scale, redshift, type of initial perturbation spectrum and statistic under consideration. For example, 
for initial spectra with significant large-scale power (see Fig. 13 in 3]) or at high redshift for CDM models, 
including just one-loop corrections to the linear power spectrum gives very good agreement with numerical 
simulations. On the other hand, for redshifts z « when the effective spectral index at the nonlinear scale 
becomes n e g ~ —1 one-loop PT deviates up to 20% for k < 0.2 /iMpc -1 (see Fig. 6 in 0). Similar results 
hold for the bispectrum 

The purpose of this paper is to cast PT in a different form, which helps to understand the physical meaning 
of different contributions in the infinite perturbation series, following ideas initially put forward in |(|. In 
particular we show that different class of diagrams organize themselves into a few characteristic quantities, 
the most important of which (in terms of improving PT) is the non linear propagator that measures the 
memory of perturbations to their initial conditions as a function of scale. When an infinite number of 
diagrams is resummed into this quantity, the remaining (infinite number of) diagrams form a well-defined 



expansion where each term with I loops corresponds to the effects of t + 1-mode-coupling, and dominates in 
a narrow range of scales. 

To illustrate these ideas and hopefully motivate the reader to go through the rest of the paper (which 
is rather technical for those not familiar with field theory methods), a quick calculation of the nonlinear 
power spectrum in the Zel'dovich approximation provides a good understanding of the difference in behavior 
at nonlinear scales between (standard) PT and the renormalized perturbation theory (RPT hereafter) we 
develop in this paper. In the Zel'dovich approximation |7J, the trajectories of particles away from their initial 
(Lagrangian) positions q are given by x = q + , 5(q) where ^(q) is the displacement field. The matter density 
perturbation, defined as S(x) = p(x)/p ~ 1, can be written in Fourier space as, 
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by means of the mass conservation relation [1 + <5(x)]<i 3 a; = d 3 q. This gives for the power spectrum, defined 
asP(Jfe)Mk + k') = (a(k)*(k / )), 
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where we introduced r = q q', A* = *(q) - ^(q') and used that (e* K '*) = 1 from Eq. (JTJ. In the 
Zel'dovich approximation the displacement satisfies — —i(k/k 2 )5L, where Sl stands for the linear density 
perturbation in Fourier space. Moreover, for Gaussian initial conditions ^ is a Gaussian random field and 
thus (e jk - A *) = e -!« k - A *) 2 >, which can be used to cast Eq. © 



as, 



P {k) = / AJL e^k-r \-[k^l-i { kX)\ _ {\ (3) 
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where /(k, r) = J d 3 q(k ■ q) 2 cos(q • r)P[ J (q)/q A and a 2 = l(k,0)/k 2 is the variance of the displacement 
field (and also the one-dimensional velocity dispersion in linear theory). Equation is the well-known 
relationship between the nonlinear and linear power spectrum in the Zel'dovich approximation Q • Although 
this provides a poor description of the power spectrum (with the nonlinear power being smaller than linear), 
having an exact solution such as Eq. © serves to illustrate the differences between PT and RPT. In PT, 
one basically expands in the amplitude of the linear spectrum, therefore 



m = / elk ' r EMh* 2 ^ - /(Mr = E (4) 

where the term with n = 1 gives the linear spectrum (tree-level, zero loops, i.e. I = 0), n = 2 gives one- loop 
{1=1) corrections, and so on. The left panel in Fig.HJshows the P^(k) for £ = 0, 1, 2, 3, for a CDM model 
with normalization a% = 0.9 and shape parameter T = 0.175, with nonlinear scale k n \ ~ 0.2 h Mpc -1 . Solid 
lines denote positive contributions, dashed lines negative ones. Note that at large enough scales only the 
linear spectrum (£ = 0) contributes, as expected, however as the nonlinear scale is approached different loop 
corrections become of the same order, with significant cancellations among them. Similar behavior has been 
seen in the case of PT for the exact dynamics, at least for some initial spectra, see e.g. 0. 

In RPT, as we discuss in detail below, one attempts (among other things) to sum an infinite class of 
diagrams for the propagator, which (as will be shown in section ITlI G 21 correspond to the factor of e~ fc "■<• 
in the Zel'dovich approximation, Eq. @- The remaining infinite series corresponds to the result 
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FIG. 1: Comparison between PT (left) and RPT (right) loop expansion for the nonlinear power spectrum in the 
Zel'dovich approximation. P n j denotes the exact result for the nonlinear power spectrum, Eq. J^J, whereas Pp^ 

(Jrpt) denotes the ^-loop correction in PT (RPT). Dashed lines denote negative values. The resummation of the 
propagator involved in RPT leads to a well-defined perturbation expansion even in the nonlinear regime, unlike PT. 



m = / it e*- t « = 1 © 

which is shown in the right panel in Fig^for ■£ = 0, 1, 2, 3. We see that as a result of this partial resumma- 
tion, the behavior of the resulting perturbation expansion is drastically altered: successive terms dominate 
at increasingly smaller scales, and there are no cancellations among them (all contributions are positive). 
Moreover, as we discuss below, RPT provides a physical understanding of the different contributions: the 
sum represents the contribution of n-modes coupling (the larger n the generation of power moves to smaller 
scales) , the decaying exponential factor represents how the amplitude and phases of the Fourier modes differ 
from linear evolution from the initial conditions. Indeed, it is easy to see that Eq. can be rewritten as 



P(k) 



fc(k-q L ..„) F„(q 1 ,...,q„) P L (qi)d 3 qi . . . P L (qn)d 3 q n , 



(6) 



where q x n = q x + . . . + q„, and F n represents the PT kernels [see Eqs. (|1011 3|l for a definition] in the 
Zel'dovich approximation [see Eq. I|38|) below]. The sum corresponds to all mode-coupling diagrams (with 
the appropriate n! weighting), with all "tadpole" diagrams in the language of standard PT [To| being 
summed to the exponential prefactor. This propagator renormalization factor is given by the square of 
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(Z><5(k) /T>5l{\s)) I D+ = exp(— fc 2 <7^/2), where V represents here a functional derivative |2J| and D + is the 
linear growth factor. The exponential decay is caused by the deviations from linear evolution (see sec- 
tion 1111 (T21 for a derivation), and the decay length (here a^ 1 ~ 0.16 /iMpc -1 ) defines a characteristic scale 
at which individual Fourier modes do not evolve as in linear PT. The negative sign contributions in PT (see 
left panel in Fig. ^) are an artifact of the power series expansion of propagator renormalization; as k increases 
these terms become large and unless one sums up the infinite series this leads to a breakdown of PT. Note 
that the n = 1 term in Eq. © gives a contribution to P(k) which is proportional to the linear spectrum, 

i.e. e~ k <t «Pl(/c), therefore if the linear spectrum were cutoff at some scale k c [-Pz,(fe) = for k > k c ] the 
nonlinear effects encoded in propagator renormalization do not lead to power generation at k > k c . It is in 
this sense that n > 2 encodes mode-mode coupling. 

In summary, the identification in RPT of physical quantities such as the propagator in the perturbation 
series helps to cure the lack of convergence of PT in the nonlinear regime. Notice also that the behavior of 
RPT in Fig. ^ is similar to the halo model description of the power spectrum, the similarities are explored 
below in section Ivl 

The rest of this paper is organized as follows. In the next section we review the standard results from 
PT and describe an alternative description necessary to develop RPT. Section II I II applies the result in 
the previous section to the calculation of statistics, in particular the power spectrum and the nonlinear 
propagator. Section IIVI derives the main results of RPT, and section [V] discusses the connection between 
RPT and the halo model. The conclusions are presented in section IVTl 

II. DYNAMICS 
A. In Standard Form 

For completeness, we review here the standard form of the equations of motion in Fourier space and their 
solutions in standard PT, before going into the new approach that helps us understand better the different 
terms in the perturbation series and thus facilitate the process of resummation. 

We study the evolution of Cold Dark Matter (CDM) in the single-stream approximation, where the 
relevant equations of motion are conservation of mass, momentum, and the Poisson equation (see e.g. [ill]). 
Consistent with this, the velocity field is assumed to be irrotational, then gravitational instability can be 
fully described in terms of the density contrast <5(x, r) = p(x,r)//3 — 1 and the peculiar velocity divergence 
6 = V-v. Defining the conformal time r = / dt/a, with a(r) the cosmological scale factor, and the conformal 
expansion rate TL = dlna/dr, the equations of motion in Fourier space become |12llla .ll4| 
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Or 



H(r)0(k, r) + ^n m (T)n 2 (T)~S(k, r) = - J d 3 hd 3 k 2 fc(k-k 1 -k 2 )/3(k 1 , k 2 ) 0(k l5 r) 0(k 2 , r). (8) 



The non linear interactions (in the right-hand side) are responsible for the coupling between different Fourier 
modes, given by 
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Equations (J7J and (JSJ are valid in an arbitrary homogeneous and isotropic universe, which evolves according 
to the Friedmann equations. For simplicity here we will restrict ourselves to an Einstein-de Sitter model for 
which £l m = 1 and Qa = 0. As a consequence, these equations can formally be solved with the following 
perturbative expansion, 



5(k,r) = ^a"(r)5 n (k), 0(k,r) = -H(r) £ a"(r)0 n (k), 



(9) 



n=l 



where on/y £/ie fastest growing mode is taken into account (below we discuss how to generalize this to the 
full time dependence) . The equations of motion, Eqs. H7I8|I determine S n (k) and 6 n (k) in terms of the linear 
fluctuations, 
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where F n and G„ are homogeneous functions of the wave vectors {q 1; . . . ,q„} with degree zero. They are 
constructed from the fundamental mode coupling functions according to the recursion relations |l3j | 



n-l 



F n (q 1 ,...,q n ) = 



, . . . , ^ 7nl 



< x (2n + 3)(n-l) 
+2 / 3(ki,k 2 )G n _ m (q m+1 ,...,qJ , 



(2n + l)a(ki,k 2 )F„_ m (q m+1 , . . . , q n ) 



(12) 



n-l 



G„( qi ,...,q„) = 



G m {<ii, ■ ■ ■ , q m ) 



3a(ki,k 2 )F„_ m (q m+1 , ...,q n ) 



where k x = q x + . . + q m , k 2 = q m+1 
the initial perturbations. 



-± (2n + 3)(n-l 
2n/3(ki,k 2 )G„_ m (q m+1 , . . . ,q„) , (13) 
-. . . + q„, k = k!+k 2 , and F 1 = G 1 = 1, with <5i(k) = 6>i(k) = S (k), 



B. In Matrix Form 



The equations of motion can be rewritten in a more symmetric form 0, Il5| by defining a two-component 
"vector" 



* a (k,ry) = (S(k, v ), -9(k,r])/n), 
where the index a — 1, 2; and we have introduced a new time variable, 



(14) 



n — In a(r). 
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corresponding to the number of e-folds of the scale factor. In the cosmology we are considering (f2 TO = 1), 
a(r) represents also the linear growth factor of perturbations in the growing mode. Using v E , a (k, z) one can 
rewrite Eqs. Q and (jSJ) together as (we henceforth use the convention that repeated Fourier arguments are 
integrated over), 



^* a (k,r ? )+n ai ,* 6 (k,j 7 )=7^ ) c (k,k 1) k 2 ) tf 6 (k 1)?7 ) * c (k 2 ,77), 



where 



-1 
-3/2 1/2 
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and "lube IS ^ ne symmetrized vertex matrix given by 
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^(kMM) = fc(k-k a -k 2 ) a(k 2 ,kx)/2, 
7^(k,k ls k a ) - foCk-ki-kaJ/JCki.ka), 



(18) 



and zero otherwise, 7^ (k,kj, kj) = 7^ 6 (k, kj, kj). An implicit integral solution to Eq. Hit) [I can be found 
by Laplace transforming in the variable 77, 



a^H * 5 (k,w) = o (k)+ 7 W(k,k 1 ,k 2 ) * 6 (k 1 ,o; 1 )* c (k 2 ,a;-a; 1 ), 



(19) 



where a (k) = x P a (k, 77 = 0) denotes the initial conditions, set when the growth factor is one, and a b (u>) = 
uj5 a b + £l a b- Multiplying by the matrix, 



1 



2uj + 1 2 
3 2ui 



(2w + 3)(w - 1) 

and performing the inversion of the Laplace transform, we finally get the formal solution 



* Q (k, rj) = 9ab(v) 06 (k) + f drf g ab (rj - rf) 7 « (k, k 1; k 2 ) * c (k 1( »/)*«i(k 3 , J?'), (20) 

Jo 

where the linear propagator g a b(v) is defined as (c > 1 to pick out the standard retarded propagator |l5j') 
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(21) 



for rj > 0, whereas g a b(v) = for 77 < due to causality, and g a b(v) ~ * &ab as 77 — » + . The linear propagator 
is the Green's function of the linearized version of Eq. I|ltj|) and describes the standard linear evolution of 
density and velocity fields from any configuration of initial perturbations. From a physical point of view, the 
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most interesting initial conditions are given by when i5(k, rj = 0) and 0(k,ri = 0) are proportional random 
fields (instead of independent), in which case we can write 



</> a (k) =u a S {k) (22) 

where u a is a two component "vector" . This covers the usual case of growing mode initial conditions for which 
u = (1,1), for which the second term in Eq. I|21|) does not contribute upon contraction of g a b with <pb (k), 
and decaying mode initial conditions for which u = (2/3, —1). For the sake of simplicity and definitcness, we 
will only contemplate the case </> a (k) = u a <$o(k) throughout this paper. 



C. Diagrammatic Representation of the Solution 



Equation (|20() for a given k-mode has a simple interpretation. The first term corresponds to the linear 
propagation from the initial conditions, whereas the second term contains information on non-linear inter- 
actions (mode-mode coupling) . This non-linear contribution comes from the interaction of all pairs of waves 
ki and k2 (whose sum is k due to translational invariance), at all intermediate times rf (with < rf < rf). 
The interaction is characterized by the vertex j^cd and then linearly evolved in time by g a b(i] — if)- As we 
will see shortly, this simple interpretation will lead, when applied recursively, to a graphical representation 
of the solution \& a in terms of diagrams that allows the resummation of the different statistical objects of 
interest. 

With the help of Eq. H2U|) . we seek for an explicit expression for ^(k, rf) in the form of a series expansion 
[see Eq. ©] 



<^(k,r0 = £*W(M), (23) 

n=l 

with [compare with Eqs. HlOlllfl ] 

*i n) (M)= / <5D(k-k 1 ...„)^ n) (k 1 ,...,k n ;, 7 )<5o(k 1 )... ( 5 (k n ), (24) 

where ki... n = ki + . . . + k n . Replacing Eqs. (|23I24|I into Eq. (|2"U)l . we find the recursion relation satisfied by 
the kernels [see Eqs. (|12I13|) ] 



^(ki,. 
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YT m =i Jo ds 9ab{v ~ s) 7^(k, ki... m , k 



T^ m h\r ■ c^l 'TV n, ~ 1Jl ') fir ■ c \ 

■J c ^l...mj b ) J d \ i ^m-\-l...m 



{n—m) , 



symmetrized 



(25) 



where the r.h.s. has to be symmetrized under interchange of any two wave vectors. For n = 1, J-^ " (n) = 
gab(v) u b- Equation 1251) reduces to the standard recursion relations given by Eqs. (|12I13|) in the limit that 
initial conditions are imposed in the infinite past, i.e. by replacing the lower limit of integration in Eq. (|25|l 
by s = -co, in which case only the fastest growing mode survives and JF(") = a n (F n ,G n ). Otherwise 
Eq. I|25|) gives the full time dependence of PT solutions, including all transients from initial conditions |15|. 

Although in principle the solution for the kernels gives us an analytic expression for ^W, it quickly 
becomes very cumbersome with increasing n, therefore in practice this is not an ideal way to proceed; this 
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is however the route taken in standard PT. Another reason apart from technical complication is the fact 
that the kernels (thought as functions of n) contain redundant information, since they all share the same 
building blocks, the vertex and linear propagator. In order to overcome these disadvantages we will introduce 
a graphical representation of the basic objects of the theory, and simple rules to combine them, such as to 
build a set of diagrams that are in one to one correspondence with each order vj/w. 

As discussed [see Eq. (|20() ]. there are only three basic building blocks, the initial field 4> a , the linear 

propagator g a b, and the vertex 7„ hc - Their graphical representation are shown in Fig. [21 Since we resort to a 
graphical representation of the basic building blocks rather than derived quantities such as the kernels (as 
in done in the diagrammatic approach in standard PT |l0lll2l[T3| '). the vertices in our diagrams are simpler 
than in standard PT, they always involve three lines (a consequence of quadratic nonlinearities) instead 
of a variable number of lines depending on the order in PT. The complication, at first sight, is that the 
diagrams carry information about time variables when the nonlinear interactions occur, and one integrates 
over these intermediate times (effectively summing over all possible interactions that happen between the 
initial conditions and the present time). This information is hidden in standard PT, where time evolution has 
been already "integrated out" (standard PT diagrams can be indeed thought as "collapsing" the diagrams 
shown in Fig. [3] below in the time direction). We will see here that this information can be taken advantage 
of, in fact making possible the process of resummation. 

<j> a (k) : -—^O g ab 0l-iV): 11 ~* J 1 ' Yabc( k ' k i k 2) = k^/ b 

a a b a ^ 

c 

k = k l + k 2 

FIG. 2: Graphical notation of the basic objects in the perturbative expansion, the initial field 4> a , the linear propagator 
g ab , and the vertex 7^.. 

To obtain the set of diagrams corresponding to the nth order term in Eq. I|23[) one has to draw all 
topologically different trees with n — 1 vertices (branchings) and n initial conditions. Each tree is obtained 
as follows: from the final time variable 77 we draw a line backwards until it reaches a vertex, where the line 
bifurcates into two branches, with each of them continuing until they reach another vertex or an initial field 
at i] = 0. If the branching is asymmetric it carries a factor of 2. This process is repeated at each vertex until 
all the branches end up in initial fields. 

Each diagram with n — 1 vertices represents an integral contributing to ty( n >. The rules to obtain these 
integrals are as follows. Each of the n initial fields <fi is characterized by a momentum k,;. Each branching 
corresponds to a vertex 7*- s )(k, ki, k^) and represents the interaction of two incoming waves ki and k2 
that couple to form one outgoing k (due to the quadratic nonlinearities in the equations of motion). Each 
interaction happens at a given time Sj (0 < Sj < rj) and conserves momentum (k = ki + k2). The branches 
correspond to linear propagators representing the linear evolution of a given mode k m from time variables Si 
until Sj. Notice that all the arrows throughout the diagram point away from the initial conditions, indicating 
the direction of forward time evolution. Finally, all the intermediate wave vectors are integrated over as well 
as all the time variables Sj, each between [0,ry]. The diagrams up to Vl/W are shown in Fig. EI Note that 
it is the existence of these simple rules for finding the nth order term, which are independent of n, that will 
allow us to considerably simplify the resummation of diagrams. 

As an example of these rules let us write explicitly the integrals corresponding to the SE^ 2 -* and 
diagrams in Fig. 
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*P (1) (k,r|) >P (2) (k,r|) ^'(k^) >P (4) (k,ri) 

FIG. 3: Diagrams up to order n — 4 in the series expansion of ^(k, rf). 



n 



= J d 3 k x J d 3 k 2 J ds g ab ( v - s) 7 «(k,k 1) k 2 ) g ce ( S )Mki) g df (s)cb f (k 2 ), (26) 
*i 3 )(k,r y ) = 2 J d 3 h J d 3 k 2 £ ds g ab ( V -s)^ s J d (kMM) 9ce(s)M*i) *f (k 2 ,s). (27) 

The diagrams we construct here are very similar to any field theory with quadratic nonlinearities , perhaps 
the closest example to ours is that of turbulence where similar methods are well known 16, 17, 13. There 
is also a recent paper []~9| that applies path-integral methods to gravitational clustering in cosmology. 

III. STATISTICS 
A. Initial Conditions 

As it stands, the integral Equation (|20|> can be thought as an equation for ^(k, 77) in the presence of an 
"external source" or forcing given by the initial conditions </> a (k) (i.e. <5o(k)), with prescribed statistics. Here 
we assume that the initial conditions are Gaussian; the statistical properties of </> a (k) are then completely 
characterized by its two-point correlator 

<<Mk) fo(k')> =<5 D (k + k') u a u b P Q {k), (28) 

where Po(k) denotes the initial power spectrum of density fluctuations. According to Wick's theorem, all 
higher order correlations of an odd number of fields vanish, whereas for an even number there are (2n — 1)!! 
contributions corresponding to all different pairings of the 2n fields, 

(^ 1 (k 1 )---^ 2 „(k 2n ))= II (<M k *) <M k i)>- (29) 

all pair associations p pairs (i,j) 

When statistics are calculated, pairs of initial fields <j) a are replaced by the initial power spectrum as one 
of the basic building blocks. It will be graphically denoted by the symbol shown in Fig. that arises from 
gluing a pair of initial fields from Fig. 
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FIG. 4: Diagrammatic notation for the initial power spectrum. 



B. Power Spectrum 

The goal is to calculate different correlation functions of the final field ^(k, rf), the simplest being its 
two-point correlator or power spectrum, 

( *„(k, i ? )* 6 (k / , rj) ) = fo(k + k')P„i,(k, rj). (30) 

Equations (|28I29() above give ensemble averages of a product of initial fields <j) a , with the result either 
vanishing or depending only on the initial power spectrum. This will allow us to develop a perturbative 
expansion for the final power spectrum in terms of integrals of the initial one, similar to the expansion for 
v P a (k, 77) in terms of ^ a (k) described in Eqs. I|23|) and 1)24(1. Replacing the series expansion of Eq. I|23|l into 

Eq. H30f). we obtain P a f,(k, 77) = X^fco -^at V s -' r l) where the f-loop term is given by (assuming Gaussian initial 
conditions) 

2£+l 

fo(k + lOPj?(M) = J2 (*i m) (k,ry)*f +2 "" l) (k',, ? )) . (31) 

m— 1 

We can now extend the diagrammatic representation of Section Til CI to describe each term in the above 
expansion. To draw each diagram that contributes to Pf$ , we put one of the tree- like diagrams for 
against one for \[/( 2£+2 ~ m ) with their initial fields facing each other. Next, we pair the initial fields in all 
possible ways and glue the pairs. Each pair of initial fields with characteristic wave numbers k and k, and 
indices c and d is then converted to an initial power spectrum fc(ki + k,-) u c UdPo(k-i). Finally, we do the 
same with all combinations of one tree of with one of +2 m ' (i.e. taking into account that after 
order n = 3 there is more than one tree at each order). Going through this process one often obtains the same 
diagram. Therefore, the final weight of a diagram is given by the numerical factor for each "independent" 
diagram due to its trees, as described in Fig. J2J, and another factor that takes into consideration the 
counting of the equal diagrams generated by the pairing process. All diagrams up to one-loop and some of 
the 29 two- loop diagrams are shown in Fig. (see [la] for a full account in the similar case of turbulence) . 

Two sets of diagrams have not been included in Fig. [3] since they vanish identically. One set is composed 
by diagrams that contain a sub-diagram linked to the main part by only one linear propagator. It can be 
shown that the sum of these fragments to all orders gives the ensemble average of the field, ( ^/ a (k, 77) ) , 
which is zero. For a similar reason the set of unconnected diagrams vanishes. 

For future use in the process of resummation, we define here the "principal cross section" of a diagram 
as the line that splits the diagram into two cutting only through initial power spectra. It represents the 
points where the two trees have been glued together. It is shown as a dashed line for the different diagrams 
in Fig. Each diagram in Fig. is in one-to-one correspondence with a certain integral. The rules for 
obtaining these integrals are the same to those described in Section III CI For example, the diagram 1 in 
Fig. |3 leads to the well known linear power spectrum for growing mode initial conditions (u a = (1, 1)) 
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FIG. 5: Diagrams for the correlation function P a b(k, 77) up to two-loops (only 7 out of 29 two-loop diagrams are 
shown here). The dashed lines represent the points at which the two trees representing perturbative solutions to ^ a 
and have been glued together. 



P o ( ° ) (k,r 7 )| dl = 5ac (r ? ) Uc P (k) u d g bd {v) — »ii(k,a) = a 2 (r)P (k), 
and the integral corresponding to diagram 2 (one loop) in the same figure gives, 



(32) 



P a J b \k,i 1 )\ d2 = 2 I ds x I ds 2 I d 3 q 3 (77- Sl )7 (s) (k,q,k-q)g( Sl )g( Sl )j a P (9) J Po(|k-q|) 

l(v - s 2 )-/ {s) (-k, q, -k + q) g(s 2 )g(s 2 ) 



(33) 



where we introduced the shorthand notations g a (s) — g a b{s)ub (= e s u a for growing mode initial conditions) 
and [57^ gg] a — 9ab%cd9c 9<i- As a further example, diagram 8 (two loops) gives, 

P Q ( b 2) (k,77)| d8 = 16 j d 3 q J d 3 pG a (k,q,p,7 7 ) G 6 (-k,-p,-q,7?) P (p) P (q) P (|k - p - q|), (34) 
where we made use of the symmetry with respect to the dashed line in Fig (J5J diagram 8, and defined, 

G a (k,q,p,77) = J dst J ds 2 g a b(vsi) 7bcd( k ' <*> k_t l) 9c(si)gdf(si-s 2 ) 7^(k-q, p, k-p-q) <? g (s 2 )<7fc(s2) 

(35) 
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C. Non Linear Propagator 



1. Definition and Diagrammatic Representation 



Nonlinearities can be thought as modifying the linear propagator defined in Eq. ((2111 , leading to propagator 
renormalization. The non-linear propagator is defined by 



Note that translation invariance requires the condition k = k', familiar from two-point statistics. Indeed, 
translation by some vector q, x — > x + q, gives G — > Gexp[i(k — k') • q]. Strictly speaking, this object is 
not a propagator in the usual sense of a Green's function, since nonlinearities are involved, i.e. it is not true 
that in the nonlinear case ^(k, 77) — G a b(k, 77)^6 (k). But G a b is the generalization of the linear propagator 
in the sense that its diagrammatic representation also has one entry and one exit arrow. It also generates 
all diagrams for the two-point function that can be divided into two subdiagrams by cutting a single initial 
power spectrum similar to what happens with g a b in the diagram for the linear power spectrum (diag. 1 
in Fig. 01 

With the help of the Eq. 1)23(1 , Eq. ((36(1 can be cast as a series expansion 



where we have explicitly separated the linear part from the non- linear contributions. The fully non- linear 
propagator represents the ensemble averaged response of the final density and velocity divergence fields 
to variations in the initial conditions. In more physical terms it quantifies how much information of the 
initial distribution of a k-mode remains in the final state at the same k, in an ensemble average sense. In 
paper II [2(|, we show that for Gaussian initial conditions (^ a (j>b) = G ac (4> c (f>b), therefore G can also be 
thought as a measure of the cross-correlation between final and initial conditions, or indeed as a genuine 
propagator (or Green's function) in two-point sense. At very large scales, we expect linear evolution to scale 
initial conditions by the growth factor; consequently, the non-linear corrections in Eq. 1(37(1 vanish as k — > 0. 
As one approaches the non-linear regime interactions gradually "erase" the initial distributions; therefore, 
we expect G(k, 77) to vanish as k — > 00. See section UlI G 21 below for the behavior of G in the Zel'dovich 
approximation, the case of the exact dynamics is presented in paper II |2p| | together with comparisons against 
numerical simulations. 

The series defining the nonlinear propagator can be represented in a diagrammatic fashion. For Gaussian 
initial conditions, only the odd terms in the expansion for ty a (k,r]) in Eq. ((37|) will contribute to G a b(k, 77). 
The nth non-linear correction corresponding to n loops is obtained by a functional derivative of vl/( 2 "+ 1 ). 
This can be done diagrammatically in a straightforward way. From each tree characterizing \I/( 2,l + 1 ) (k, 77) 
in Fig. |3 there will be In + 1 contributions, each obtained after dropping one initial condition represented 
by a blank circle (this accounts for the functional derivative) . The linear propagator left at that branch will 
carry a wave number k. Finally the ensemble average of the remaining 2n 0's is done by pairing the fields 
according to Gaussianity. 

Figure shows all diagrams up to two loops for the non-linear propagator. Notice that each diagram 
has one entering arrow and one exiting, connected by a chain of linear propagators that runs through the 
diagram without intersecting initial conditions. This path is called the 11 principal path" , and will be used 
later in the process of resummation. 




(36) 




(37) 



12 




FIG. 6: Diagrams for the propagator G(k, if) up to two loops. 



2. A simple example: Zel'dovich Approximation 

To illustrate the behavior of the nonlinear propagator we go back to the discussion in the introduction 
about the Zel'dovich Approximation. Let us work out the density propagator, G$ = Gn + G22 which 
assuming growing mode initial conditions follows directly from taking the derivative of the final density with 
respect to the initial one. For simplicity we consider only the fastest growing mode at each order, so we can 
use the standard PT kernel representation, Eqs. JJjJ, 110(1 and ((12(1 . and take advantage of the simple form 
the density field kernel takes in this approximation |2l| 

F n (q 1 ,... ! qJ = i^|li...^ ) (38) 

where k = q x + . . . + q n . Due to Gaussian initial conditions, only odd terms in the series expansion for 8 
contribute, therefore 



G s = a ^(2n+l)!! / F 2n+ i(k, q x , -q 1; . . . , q„, -q n )P L (qi, a)d 3 qi . . . P L (q n , a)d 3 q n , (39) 

n=0 

where the term n — gives just the scale (and growth) factor. Using Eq. 1(38(1 this series can be summed up 
right away, 



Gs = a 



^ nil 2 / 



k-q 

2 



P L (q,a)d i q = a exp(-fcV^/2), 



(40) 



where <r\ = (1/3) J d 3 qPL(q,a)/q 2 . The extension of this result to the exact dynamics for density and 
velocity fields and comparison with measurements in numerical simulations is given in pfij . 
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D. Vertex 



The vertex defined in Eq. (|18f) is the third basic object of the theory, and is also renormalized due to 
higher-order nonlinear interactions. We define the symmetric full vertex through the relation 



<*y,(M) 

<ty e (ki)<ty/(k 2 ) 



P7] PS PS 

2 ds dsi ds 2 G ab (t] - s) Iy d (k, s; ki, si ; k 2 , s 2 ) G ce (si) G df (s 2 ), (41) 
Jo Jo Jo 



thus r^(ks,k 1 s 1) k 2 S2) = 7bcd(k, k x , k 2 ) <5 D (s - s x ) <5 D (s - s 2 ) + "higher-order corrections". 

The full vertex depends not only on the two incoming wave vectors but also on the corresponding times 
when they "enter" the interaction. It can be shown that the outgoing wave vector k is determined by wave 
vector conservation, r^(k, s; ki, si; k 2 , s 2 ) oc <5o(k — ki — k 2 ). The diagrammatic representation up to one 
loop is shown in Fig. 




+ 




+ 




FIG. 7: Diagrams corresponding to the renormalized vertex F, Eq. 1411 . up to one loop. 



IV. RENORMALIZED PERTURBATION THEORY (RPT) 

A. Non-Linear Propagator 

In this section we will perform the first step in the resummation of the series in Eq. (O. We divide the 
diagrams in Fig. into two classes. The Principal Path Reducibles are the ones that can be split in two 
disjoint pieces by cutting one linear propagator belonging to the principal path. Diagram 3 in Fig. is 
the only example up to two loops. Note that the linking propagator represents the linear evolution of the 
same Fourier mode that enters and leaves the diagram, and hence it corresponds to the "independent" wave 
vector that is not integrated over. Diagrams that cannot be split in this way are denoted as Principal Path 
Irreducibles. 

All the Principal Path Irreducible diagrams start with a linear propagator g(s 2 ) entering the diagram at 
s 2 and finish with g(n — s±) exiting at si, with an irreducible structure in between. We define as S(k, s±, s 2 ) 
the sum of all these irreducible structures. Thus, we have 



Sum of all irreducibles = g a b{v) + / / ds 

Jo Jo 



2 gaciv - si)S cd (k, sx,s 2 ) gdb(s 2 ) 



(42) 



The relation r\ > si > s 2 > follows from causality along the principal path. All the diagrams for 
S(k, s%, s 2 ) up to two loops are shown in Fig. EI 
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FIG. 8: Diagrams for S(k, 81,82) up to two loops. The dashed lines enclose reducible fragments, see section llV CI 



The Reducible diagrams are built by two independent fragments linked together by one linear propagator 
that evolves up to some time, say s 2 . The diagrams also finish with a linear propagator g(rj — s\). This can 
be clearly seen in diagram 3 in Fig. There we see that between si and s 2 there is the first contribution to 
S(k, si, S2) in Fig. El Before s 2 , there is the first non-linear term in the expansion for G(k, s 2 ), that is diagram 
2 in the same Figure. The higher-order Reducible diagrams will have all the higher-order contributions to £ 
between s\ and s 2 , and before S2 all the terms in the series for G except the linear one, g. We then write, 

rv r»i 

Sum of all reducibles = / ds\ / ds 2 g ac (v - s%) £ C( j(k, Si, s 2 ) [Gdb(k, s 2 ) - gdb(s 2 )} (43) 
Jo Jo 

Adding up Eqs. (|42l) and 143fl . we arrive at the equation for the non- linear propagator (usually known as 
Dyson's equation) 

rn rsi 

G ab (k,r)) = gabWl + / ds x / ds 2 g ac {r] ~ s\) S cc j(k, si,s 2 ) G d6 (k, s 2 ). (44) 
Jo Jo 

For further reference we generalize the previous equation to the case where the non-linear propagator 
evolves from a general time variable r( ^ 0, 

G a b(k,?7, 77') = g afc (7? - 7/) + / dsi / ds 2 g ac (?7- si)E cd (k, si,s 2 )G db (k, s 2 ,7/). (45) 

J 7/' J 77' 

Notice that G fc(k, 77) = G &(k, 77, 77' = 0). We will represent the non-linear propagator G by a thick 
straight line, shown in Fig. I1UI together with the graphical representation of Eq. (|45|) . 
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B. Power Spectrum 



We now turn to the derivation of the equivalent of Eq. (|44() for the non-linear power spectrum. Let us 
first consider in Fig. [S] all the diagrams that contain only one initial power spectrum Po(^) at the "principal 
cross section". Examples of this type are diagrams 1,3,4 and 7. They all share the property that at each 
side of the "principal cross section" there is an independent sub-diagram with one entry and one exit arrow 
connected by a "principal path" . By summing diagrams of this sort to all orders one ensures to account for 
all the contributions to the non-linear propagator at both sides of Po(k) at the "principal cross section" 

Sum of all diagrams with one initial power 

spectrum at the "principal cross section" = G ac (k, rf) u c Po(k) Ud Gbd(— k, if). (46) 

Consider next diagrams containing two or more initial power spectra at their "principal cross section" . 
From those, diagrams 2,8,9,10,11 in Fig. [5] form a special class since they cannot be split in two parts 
(with more than one g each) by cutting one linear propagator. They are all characterized by an irreducible 
structure enclosing the "principal cross section" finishing with two vertices at, say, S\ to the left and s 2 to 
the right. They continue after each vertex with a linear evolution from Si to r\ with g{r\ — si). We define 
as $(k, si, S2) the sum of all these irreducible structures. All diagrams for $ up to two loops are shown in 
Fig. 

The remaining diagrams also contain a contribution to $(k, s\,S2) around the "principal cross section", 
but either to the left of si or to the right of S2 they have a non-linear contribution to the non-linear propagator 
G(k, 77, Si) (diagrams 5 and 6 in Fig. [3] respectively). Contributions to the non-linear propagator at both 
sides of <£> appear for the first time in three-loop diagrams. 

After adding the diagrams described above to all orders, we arrive at 

Sum of all diagrams with two or more initial power (47) 

/•J7 rri 

spectra at the "principal cross section" = / dsi / ds2 G ac (k, 77, si) $ c d(k, si, s 2 ) G b d(— k, rj, s 2 ). 

Jo Jo 

Adding Eqs. (|46I47|) we arrive at the integral equation for the non-linear power spectrum, 

rv rn 

P a b(k,r)) = G ac Qs.,r))u c P (k)udG bd (-'k,r])+ / dsi / ds 2 G ac (k, r], si) $ crf (k, si, s 2 ) G bd (-k, ?/, s 2 ). (48) 

Jo Jo 

For further reference we extend the previous analysis to the case of a two-point correlation of final fields 
at different times 

(iMM) ^&(k',7/)> = <fo(k + k') P ab (k,T],r)'). (49) 
The equivalent of Eq. (|48|l for P a b(k,r],r]') reads 

rV rn' 

Pab(k,r/,r)') = Gac(Kil) u cPo(k)udG bd (-^il') + / dsi / ds 2 G ac (k, 77, si) $ crf (k, si, s 2 ) G bd (-k, r/, s 2 ). 

Jo Jo 

(50) 
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FIG. 9: All diagrams for $ up to two loops. The dashed lines enclose reducible fragments for propagator and power 
spectrum renormalization. The dash-dotted lines enclose reducible fragments that lead to vertex renormalization, see 
section llV (Jl for a detailed discussion. 



G ab (k,r|,r|') 



rp^! 1 (k ' S l' S 2) Ip%< 



P, h Gwl') 



- m ^<E>(k,Sj,s 2 ) -, 



ri n n n 

FIG. 10: Integral Equations 145H and 1501 for the non-linear propagator and the non-linear power spectrum. 



This equation is shown diagrammatically in Fig. I1UI where we represent the non-linear power spectrum 
by a thick line with a filled circle at the middle. 

It is beyond of the scope of this paper to consider higher-order correlations in detail, but Fig. ^] shows 
the one-loop diagrams (after renormalization) for the bispectrum 

(* Q (ki,77)* 6 (k 2 ,7?)*c(k3,?7)> = fo(k 1 +k a +k 3 ) B abc (k u k 2 , k 3 , 77), (51) 
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C. Further Resummations 



In Sections IIV Al and IIVBI we derived equations that allow the renormalization of the propagator and 
power spectrum given by the kernels £ and <F However the kernels themselves were functions of the 
linear propagator and linear power spectrum. In this section we complete the resummation procedure by 
renormalizing £ and $, which means that we find expansions for these quantities in terms of non- linear 
propagators and non- linear power spectra. We will also carry out a vertex renormalization for <E>. 

We begin by noticing that all diagrams for both £ and $, in Figs.|SJand|!5|respectively, contain two ending 
nodes. For £ there is one entry point and one exit, while for $ both are exits. In those figures we define 
as a "reducible fragment" any part of a diagram that can become disjoint from the main part by cutting 
two linear propagators. We refer as "main part" of a diagram to the one containing the two ending points 
and consider the two linking propagators as belonging to the reducible fragment. Reducibles fragments are 
shown, enclosed by dashed lines, in diagrams 2,3,4,5 for £ in Fig. [5J and in diagrams 2, 3, 4 for $ in Fig. [5] 

Notice that there are two types of reducibles fragments, which we call type G and type P, since they lead 
to propagator renormalization (G) and power spectrum renormalization (P). Type-G fragments contain a 
principal path, with one linear propagator entering the fragment and another exiting it. Examples are the 
fragments in diagrams 2,3,4 for £ and 3,4 for $. Type-P fragments contain a local principal cross section, 
with two linear propagators going out of the fragment, like the one in diagrams 5 for £ and 2 for 

Let us analyze first diagrams with type-G fragments. In the series for £, we see that diagram 2 can be 
generated by replacing g(si — S2) in diagram 1 by the first non-linear contribution to G(k, si, S2) (shown 
in diagram 2 in Fig. |SJ). Likewise diagrams 3,4 can be obtained out of diagram 1 by replacing linear by 
non-linear propagators. In the series for $, diagrams 3 and 4 can be obtained by substituting one linear 
propagator in diagram 1 by the first non-linear correction to G(k, s). 

Next, we turn to diagrams containing type P fragments. In Fig. [SJ diagram 5 can be obtained by replacing 
the linear power spectrum in diagram 1 by the first non-linear contribution to P(k, s±, S2) (shown in diagram 
2 of Fig. UJl. Similarly, in Fig. El diagram 2 is equivalent to diagram 1 with the linear propagator replaced by 
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the one-loop correction. 

It is possible to show that the sum of type-G fragments is the complete non-linear correction to the 
non-linear propagator, while the sum of type-P fragments equals the complete non-linear correction to 
the non-linear power spectrum. This completes the resummation procedure for E in terms of non-linear 
propagators and power spectra. The resummed diagrams, up to two loops, are shown in Fig. 

The renormalization of $ can be taken one step further by performing the resummation of the vertex. 
Diagrams 6 and 7 in Fig. [21 can be obtained from diagram 1 by replacing the corresponding tree- level vertex 
7 by the first one-loop contribution to the full vertex, T, shown in Fig. [7J Similarly, diagrams 8 and 9 are 
reproduced with the help of the second and third one-loop diagrams in Fig. [7J 

With the renormalization procedure described in this paper for power spectrum, propagator and vertex 
we were able to account for all the diagrams in the series for <J> up to two loops, with the renormalized series 
given in Fig ll3l In [16(, the author carried out the same procedure one order higher, and could describe all 
the series for $ up to three loops, with 19 new three-loop irreducible diagrams added to those in Fig. IT51 
However, in |l7lll8| it has been pointed out the necessity to introduce three different vertices to account for 
all diagrams, in the series of both £ and after two loops. A detailed discussion of vertex renormalization 
is beyond the scope of this paper, and will be discussed elsewhere. In our companion paper, [2(j, we get 
around this issue for the case of the non linear propagator (i.e. S) by summing the bare perturbation series 
in Fig. EJ 

In summary, it is possible to obtain the series for $, at least up to two loops, with only the corresponding 
irreducibles diagrams written in terms of non-linear propagators, power spectra, and full vertices. The first 
diagram in Fig.|5]is the irreducible at the one-loop level, while the fifth is the only irreducible at the two- loop 
level. Figure IT51 shows the final "renormalized" results for <I>. 

It is worth mentioning that the same resummation procedure can be done for the full vertex in Fig. [7J The 
"renormalized" series for is, to one loop, still given by the diagrams in Fig. [7J but with all quantities 
changed by their fully non-linear counterparts pl| . 



The expression in Eq. 148(1 corresponding to the bottom diagram in Fig.^Jhas some interesting similarities 
with the expression for the power spectrum in the halo model 22]. In RPT, the power spectrum is written 
as the sum of two terms, one that dominates at large scales (similar to the 2-halo term) and the other 




+ 




FIG. 12: Renormalized series expansion for E up to two loops. 



V. 



RPT AND THE HALO MODEL 
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FIG. 13: Renormalized series expansion for $ up to two loops. 



(involving <!>) that dominates at small scales (similar to the 1-halo term). The RPT expression for the 
Zcl'dovich approximation, shown in the right panel in Fig. ^ and Eq. © makes the analogy very clear: the 
n = 1 term corresponds to the linear power spectrum modulated by a scale-dependent "bias factor" due 
to the renormalized propagator that becomes suppressed in the nonlinear regime; whereas the n > 2 sum 
(weighted by the square of the PT kernels) resembles the sum over halo masses (weighted by the square of 
the halo profiles) in the 1-halo power spectrum. 

The analogy can be taken a bit further by calculating the nonlinear propagator in the halo model, where the 
density can be written as a sum over halos of mass to, and profile u mi (x) (normalized so that J d 3 xu m (x) = 1 ) 



P(x) 



E 



mi u mi (x - 



J m dm d 3 x' u m (x — x') S^(m — m i) 5d( x i — x ')> 



(52) 



where the mass function is obtained from the average Q^. <5d(ti — m.i)S-£,(xi — x')} = n(m) and the mean 
density is p — J m dm n(m). As we show in [2(J the nonlinear propagator can also be calculated for Gaussian 
initial conditions from cross-correlating final and initial fields, and for growing mode initial conditions the 
density propagator becomes Gg(k) = (<5(k)<$o(k'))/(<5o(k)5o(k')}, which reduces to the growth factor D + in 
linear theory (Sl — D + 5 Q ). The correlation between halo positions and the initial density field is not known 
in detail, but following the standard assumptions in the calculation of the 2-halo term, we can model it 
introducing the halo linear bias factor b\(m) 



<fo(m - mi) <5 D (xi - x') ~ n(m) 1 + 6i(m) <5l(x') 

i 

which implies, after straightforward algebra, that in the halo model, 



(53) 



D f 

Gs(k) = J mdmn(m)bi(m) u m (k) = D±b%(k), (54) 

in other words, propagator renormalization in RPT exactly corresponds to the standard (scale-dependent) 
bias in the framework of the halo model, making the analogy between the 2-halo term and the first term in 
Eq. H48|) exact. In practice the suppression from the propagator in RPT is stronger than that in the halo 
model [due to the halo profile in Eq. l|54Jl ] , however this is a simplification of Eq. H53J) where nonlinear halo 
bias and exclusion effects are ignored. 

Regarding the 1-halo term, it is difficult to connect it formally to the mode-coupling kernel $ in Eq. (|48|) 
because there is no perturbative description of such a term and $ does not appear to have a simple expression 
apart from obeying Eq. I)48[l. However, physically, the relationship is obvious: the second term in Eq. (|48|l 
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describes mode-mode coupling and so does the 1-halo term (given the usual assumptions where the 2-halo 
term is proportional to the linear spectrum). In this sense the correspondence between the description of 
the power spectrum in RPT and the halo model is exact, although the ingredients are very different. 

VI. CONCLUSIONS 

We developed a new way of looking at cosmological perturbation theory, which makes possible a well- 
controlled description of gravitational clustering at non linear scales. This formalism follows the growth 
of structure as it develops in time, decomposed into linear propagation plus interactions, summing over 
all possibilities. Although at first sight such time-decomposition seems more complicated than in standard 
perturbation theory (where time evolution is already integrated out up to the present time), we showed that 
one can take advantage of this extra information by finding (based on topological considerations) infinite 
subset of diagrams that can be resummed and identified with physical quantities. 

The most important aspect of this resummed theory, which we call renormalized perturbation theory 
(RPT), is that the nonlinear (renormalized) propagator (which enters as a key ingredient in the calculation 
of correlation functions) is a strong function of scale, decaying nearly exponentially at nonlinear scales, 
see |2(j for measurements in numerical simulations. When the perturbative expansion of the nonlinear 
power spectrum is divided into terms that sum up to the renormalized propagator plus terms that describe 
mode-mode coupling, the resulting (renormalized) perturbation theory is well behaved in the sense that each 
term dominates in a narrow range of scales and is suppressed otherwise. We also showed that this description 
turns out to be analogous to that in the halo model. 

Future work is needed to decide whether this approach will be fruitful in quantitative terms; for example, 
by comparing predictions for the nonlinear power spectrum against N-body simulations. In a companion 
paper |20| we present a detailed analysis of the resummation of the propagator and comparison of it against 
numerical simulations. The calculation of the nonlinear power spectrum in RPT requires the inclusion of 
many terms (loops) in the mode-coupling series to cover an interesting range of scales, which at present 
appears as a nontrivial task. However, there are symmetries (e.g. Galilean invariance, see ^3) that connect 
the resummation of the mode-coupling series with that of the propagator, which one might be able to take 
advantage of. This issue deserves further work and will be discussed elsewhere [23|. 

Finally, one might wonder whether the single-stream approximation made so far in our formalism will 
break down before one can get any interesting results into the nonlinear regime |25|. Although this is a 
possibility, an educated guess is that such effects do not show up until scales well in the nonlinear regime, if 
one uses as a proxy for when multistreaming becomes important the scale at which the power in the vorticity 
of the velocity field equals that in the divergence . This happens (not surprisingly, since vorticity should 
slow down the growth of power) at roughly the same scale of the "virial turnover" of the nonlinear power 
spectrum, i.e. k ~ 1 /iMpc" 1 for CDM models at z — 0. Having a robust prediction of the nonlinear power 
spectrum as a function of cosmological parameters up to such scales would be very useful for applications in 
many aspects of large-scale structure, e.g. weak gravitational lensing surveys. 
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